# B1
source('scripts/createDataset.R')

slope_temp <- 'separate'
order_temp <- 1

dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.))
reg1 <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp)

contVar <- cont1
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg2 <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, 
                                 covariates = contVar %>% paste(., collapse = '+'))

contVar <- cont2
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg3 <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, 
                                 covariates = contVar %>% paste(., collapse = '+'))

contVar <- cont3
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg4 <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, 
                                 covariates = contVar %>% paste(., collapse = '+'))

contVar <- cont4
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg5 <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, 
                                 covariates = contVar %>% paste(., collapse = '+'))

contVar <- cont5
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg6 <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, 
                                 covariates = contVar %>% paste(., collapse = '+'))

contVar <- cont6
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg7 <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, 
                                 covariates = contVar %>% paste(., collapse = '+'))

contVar <- cont7
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg8 <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, 
                                 covariates = contVar %>% paste(., collapse = '+'))

# STORE RESULTS
fileName <- paste0('results/APPENDIX_B/B1.html')
reg <- list(reg1, reg2, reg3, reg4, reg5, reg6, reg7, reg8)
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
modelsummary(reg, vcov = robust_temp, exponentiate = F, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, 
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

############

# B2
source('scripts/createDataset.R')

contVar <- contVarBase
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))

reg1 <- dat_temp1 %>% rdd_reg_lm(slope = 'same', bw = bw, order = 1, covariates = contVar %>% paste(., collapse = '+'))
reg2 <- dat_temp1 %>% rdd_reg_np(slope = 'separate', bw = bw, covariates = contVar %>% paste(., collapse = '+'))
reg3 <- dat_temp1 %>% rdd_reg_np(slope = 'same', bw = bw, covariates = contVar %>% paste(., collapse = '+'))
reg4 <- dat_temp1 %>% rdd_reg_lm(slope = 'separate', bw = bw, order = 2, covariates = contVar %>% paste(., collapse = '+'))
reg5 <- dat_temp1 %>% rdd_reg_lm(slope = 'same', bw = bw, order = 2, covariates = contVar %>% paste(., collapse = '+'))
reg6 <- dat_temp1 %>% rdd_gen_reg(slope = 'separate', bw = bw, order = 1, fun = glm, family = binomial(link = 'logit'),
                                  covariates = contVar %>% paste(., collapse = '+'))
reg7 <- dat_temp1 %>% rdd_gen_reg(slope = 'same', bw = bw, order = 2, fun = glm, family = binomial(link = 'logit'),
                                  covariates = contVar %>% paste(., collapse = '+'))
reg8 <- dat_temp1 %>% rdd_gen_reg(slope = 'separate', bw = bw, order = 2, fun = glm, family = binomial(link = 'logit'),
                                  covariates = contVar %>% paste(., collapse = '+'))

# OLS
fileName <- paste0('results/APPENDIX_B/B2_1.html')
reg <- list(reg1, reg4, reg5)
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
modelsummary(reg, vcov = robust_temp, exponentiate = F, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, 
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

# NON-PARAMETRIC CANNOT BE EXPORTED

# LOGIT
fileName <- paste0('results/APPENDIX_B/B2_2.html')
reg <- list(reg6, reg7, reg8)
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
modelsummary(reg, vcov = robust_temp, exponentiate = T,
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, 
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

############

# B3
source('scripts/createDataset.R')

dat_temp1 <- rdd_data(df$y, df$x, covar = df[,c('pid', contVar)], cutpoint = cutoff) %>% subset(., complete.cases(.))
reg1 <- dat_temp1 %>% rdd_reg_lm(slope = 'same', bw = bw, order = 1, 
                                 covariates = contVar %>% paste(., collapse = '+'))
reg2 <- dat_temp1 %>% rdd_reg_lm(slope = 'separate', bw = bw, order = 1, 
                                 covariates = contVar %>% paste(., collapse = '+'))
reg3 <- dat_temp1 %>% rdd_reg_lm(slope = 'same', bw = bw, order = 2, 
                                 covariates = contVar %>% paste(., collapse = '+'))
reg4 <- dat_temp1 %>% rdd_reg_lm(slope = 'separate', bw = bw, order = 2, 
                                 covariates = contVar %>% paste(., collapse = '+'))

fileName <- paste0('results/APPENDIX_B/B3.html')
reg <- list(reg1, reg2, reg3, reg4) %>% lapply(., function(x) as.lm)
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
reg <- list(reg1, reg2, reg3, reg4)
modelsummary(reg, 
             vcov = list(vcovCL_custom(reg1), vcovCL_custom(reg2), vcovCL_custom(reg3), vcovCL_custom(reg4)), 
             exponentiate = F, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, 
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

############

# B5
source('scripts/createDataset.R')

slope_temp <- "separate"
order_temp <- 1

dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.)) %>% mutate(D = as.numeric(x >= 0))
reg1 <- lm(formula = 'y~D', data = dat_temp1)

contVar <- cont1
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.)) %>% mutate(D = as.numeric(x >= 0))
reg2 <- lm(formula = paste0('y~D+', paste(contVar, collapse = '+')), data = dat_temp1)

contVar <- cont2
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.)) %>% mutate(D = as.numeric(x >= 0))
reg3 <- lm(formula = paste0('y~D+', paste(contVar, collapse = '+')), data = dat_temp1)

contVar <- cont3
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.)) %>% mutate(D = as.numeric(x >= 0))
reg4 <- lm(formula = paste0('y~D+', paste(contVar, collapse = '+')), data = dat_temp1)

contVar <- cont4
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.)) %>% mutate(D = as.numeric(x >= 0))
reg5 <- lm(formula = paste0('y~D+', paste(contVar, collapse = '+')), data = dat_temp1)

contVar <- cont5
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.)) %>% mutate(D = as.numeric(x >= 0))
reg6 <- lm(formula = paste0('y~D+', paste(contVar, collapse = '+')), data = dat_temp1)

contVar <- cont6
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.)) %>% mutate(D = as.numeric(x >= 0))
reg7 <- lm(formula = paste0('y~D+', paste(contVar, collapse = '+')), data = dat_temp1)

contVar <- cont7
dat_temp1 <- rdd_data(df$y, df$x, covar = df[,c('pid', contVar)], cutpoint = cutoff) %>% subset(., complete.cases(.)) %>% mutate(D = as.numeric(x >= 0))
reg8 <- lm(formula = paste0('y~D+', paste(contVar, collapse = '+')), data = dat_temp1)

fileName <- paste0('results/APPENDIX_B/B5.html')
reg <- list(reg1, reg2, reg3, reg4, reg5, reg6, reg7, reg8) %>% lapply(., function(x) as.lm)
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
reg <- list(reg1, reg2, reg3, reg4, reg5, reg6, reg7, reg8)
modelsummary(reg, 
             vcov = robust_temp, 
             exponentiate = F, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, 
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

############

# B6-B8
source('scripts/createDataset.R')

slope_temp <- "separate"
order_temp <- 1

# LINEAR
month_seq <- seq(60, 240, 12)
reg <- list()
for(b in 1:length(month_seq)){
  dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVarBase], cutpoint = cutoff) %>% subset(., complete.cases(.))
  reg[[b]] <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = month_seq[b], order = order_temp, 
                                       covariates = contVarBase %>% paste(., collapse = '+'))
  reg[[b]]$param <- month_seq[b]
}

fileName <- paste0('results/APPENDIX_B/B7&8.html')
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
modelsummary(reg, vcov = robust_temp, exponentiate = F, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, 
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

p1 <- reg %>% 
  extractRegResultsForPlot(reg = ., robust_temp = robust_temp, exponentiate = F) %>%
  createATTPlotBW(., ybreaks = seq(-2.5, 2.5, 0.05), ylimits = c(0, 0.25), yintercept = 0)
fileName_plot <- paste0('results/APPENDIX_B/B6_A.png')
ggsave(fileName_plot, p1, width = 7, height = 2.5)

# NONPARAMETRIC
month_seq <- seq(60, 240, 12)
reg <- list()
for(b in 1:length(month_seq)){
  dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVarBase], cutpoint = cutoff) %>% subset(., complete.cases(.))
  reg[[b]] <- dat_temp1 %>% rdd_reg_np(slope = slope_temp, bw = month_seq[b], 
                                       covariates = contVarBase %>% paste(., collapse = '+'))
  reg[[b]] <- reg[[b]][['RDDslot']][['model']]
  reg[[b]]$param <- month_seq[b]
}

p1 <- reg %>% 
  extractRegResultsForPlot(reg = ., robust_temp = robust_temp, exponentiate = F) %>%
  createATTPlotBW(., ybreaks = seq(-2.5, 2.5, 0.05), ylimits = c(0, 0.25), yintercept = 0)
fileName_plot <- paste0('results/APPENDIX_B/B6_B.png')
ggsave(fileName_plot, p1, width = 7, height = 2.5)

# LOGIT
month_seq <- seq(60, 240, 12)
reg <- list()
for(b in 1:length(month_seq)){
  dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVarBase], cutpoint = cutoff) %>% subset(., complete.cases(.))
  reg[[b]] <- dat_temp1 %>% rdd_gen_reg(slope = slope_temp, bw = month_seq[b], order = order_temp, 
                                        fun = glm, family = binomial(link = 'logit'),
                                        covariates = contVarBase %>% paste(., collapse = '+'))
  reg[[b]]$param <- month_seq[b]
}

p1 <- reg %>% 
  extractRegResultsForPlot(reg = ., robust_temp = robust_temp, exponentiate = T) %>%
  createATTPlotBW(., ybreaks = seq(-5.5, 5.5, 0.25), ylimits = c(0.9, 3), yintercept = 1)
fileName_plot <- paste0('results/APPENDIX_B/B6_C.png')
ggsave(fileName_plot, p1, width = 7, height = 2.5)

############

# B9
source('scripts/createDataset.R')

parties <- c("partyGRUNE", "partySPD", "partyFDP", "partyLINKE", 'partyAFD')
slopevec <- "separate"
ordervec <- 1
reg <- list()
df1 <- df

for(p in 1:length(parties)){
      
  party_temp <- parties[p]
  df <- df1 %>% select(-y) %>% rename(y = party_temp)

  dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
  reg[[p]] <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, 
                                       covariates = contVar %>% paste(., collapse = '+'))
  
}

fileName <- paste0('results/APPENDIX_B/B9.html')
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
modelsummary(reg, vcov = robust_temp, exponentiate = F, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, 
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

############

# B10
source('scripts/createDataset.R')

parties <- c('HHpartyUNION', "HHpartyGRUNE", "HHpartySPD", "HHpartyFDP", "HHpartyLINKE", 'HHpartyAFD')
slopevec <- "separate"
ordervec <- 1
reg <- list()
df1 <- df

for(p in 1:length(parties)){
  
  party_temp <- parties[p]
  df <- df1 %>% select(-y) %>% rename(y = party_temp)
  
  dat_temp1 <- rdd_data(df$y, df$x, covar = df[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
  reg[[p]] <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp, 
                                       covariates = contVar %>% paste(., collapse = '+'))
  
}

fileName <- paste0('results/APPENDIX_B/B10.html')
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
modelsummary(reg, vcov = robust_temp, exponentiate = F, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, 
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

############

# B11
source('scripts/createDataset.R')

mediatorvec <- c("age_olderthan47", "nchild_lessthan2", "nchild_morethan1", 
                 "married_married", "married_notmarried", "income2_lessthan850")
slope_temp <- "separate"
order_temp <- 1

reg <- list()
dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.))
reg_baseline <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp)
for(m in 1:length(mediatorvec)){
  mediator_temp <- mediatorvec[m]
  dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.))
  df1 <- df %>% subset(., get(mediator_temp) == 1)
  dat_temp1 <- rdd_data(df1$y, df1$x, covar = df1[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
  reg[[m]] <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp)
  reg[[m]]$pval <- t.test2(reg_model1 = reg_baseline, reg_model2 = reg[[m]]) %>% .[['p.value']] %>% round(., 4)
  reg[[m]]$baseline <- mediator_temp
}
fileName <- paste0('results/APPENDIX_B/B11.html')
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
info <- data.frame() %>% 
  rbind(c('p-value', unlist(lapply(reg, function(x) x$pval)) %>% as.character())) %>%
  rbind(c('sub-sample', unlist(lapply(reg, function(x) x$baseline)) %>% as.character()))
modelsummary(reg, vcov = 'HC1', exponentiate = F, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, add_rows = info,
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

############

# B12
source('scripts/createDataset.R')

mediatorvec <- c("age_olderthan47", "nchild_lessthan2", "nchild_morethan1", 
                 "married_married", "married_notmarried", "income2_lessthan850")
slope_temp <- "separate"
order_temp <- 1

reg <- list()
dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.))
reg_baseline <- dat_temp1 %>% rdd_gen_reg(slope = slope_temp, bw = bw, order = order_temp,
                                          fun = glm, family = binomial(link = 'logit'))
for(m in 1:length(mediatorvec)){
  mediator_temp <- mediatorvec[m]
  dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.))
  df1 <- df %>% subset(., get(mediator_temp) == 1)
  dat_temp1 <- rdd_data(df1$y, df1$x, covar = df1[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
  reg[[m]] <- dat_temp1 %>% rdd_gen_reg(slope = slope_temp, bw = bw, order = order_temp,
                                        fun = glm, family = binomial(link = 'logit'))
  reg[[m]]$pval <- t.test2(reg_model1 = reg_baseline, reg_model2 = reg[[m]]) %>% .[['p.value']] %>% round(., 4)
  reg[[m]]$baseline <- mediator_temp
}
fileName <- paste0('results/APPENDIX_B/B12.html')
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
info <- data.frame() %>% 
  rbind(c('p-value', unlist(lapply(reg, function(x) x$pval)) %>% as.character())) %>%
  rbind(c('sub-sample', unlist(lapply(reg, function(x) x$baseline)) %>% as.character()))
modelsummary(reg, vcov = 'HC1', exponentiate = T, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, add_rows = info,
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

############

# B13
source('scripts/createDataset.R')
df <- df %>% 
  mutate(turnout09_cs = replace_na(turnout09_cs, 0),
         turnout13_cs = replace_na(turnout13_cs, 0),
         turnout17_cs = replace_na(turnout17_cs, 0)) %>%
  mutate(turnout09_yes = as.numeric(turnout09_cs == 1),
         turnout13_yes = as.numeric(turnout13_cs == 1),
         turnout17_yes = as.numeric(turnout17_cs == 1),
         turnout09_no = as.numeric(turnout09_cs == 0),
         turnout13_no = as.numeric(turnout13_cs == 0),
         turnout17_no = as.numeric(turnout17_cs == 0)) %>%
  mutate(politmember_yes = as.numeric(politmember == 1),
         politmember_no = as.numeric(politmember != 1))

mediatorvec <- c('turnout09_yes', 'turnout09_no', 'turnout13_yes', 'turnout13_no', 
                 'politmember_yes', 'politmember_no')
slope_temp <- "separate"
order_temp <- 1

reg <- list()
dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.))
reg_baseline <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp)
for(m in 1:length(mediatorvec)){
  mediator_temp <- mediatorvec[m]
  dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.))
  df1 <- df %>% subset(., get(mediator_temp) == 1)
  dat_temp1 <- rdd_data(df1$y, df1$x, covar = df1[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
  reg[[m]] <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp)
  reg[[m]]$pval <- t.test2(reg_model1 = reg_baseline, reg_model2 = reg[[m]]) %>% .[['p.value']] %>% round(., 4)
  reg[[m]]$baseline <- mediator_temp
}
fileName <- paste0('results/APPENDIX_B/B13.html')
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
info <- data.frame() %>% 
  rbind(c('p-value', unlist(lapply(reg, function(x) x$pval)) %>% as.character())) %>%
  rbind(c('sub-sample', unlist(lapply(reg, function(x) x$baseline)) %>% as.character()))
modelsummary(reg, vcov = 'HC1', exponentiate = F, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, add_rows = info,
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

############

# B14
source('scripts/createDataset.R')
pid_swing1 <- df %>% subset(., !is.na(lrscale) & syear %in% c(2014) & lrscale %in% c(4:6)) %>% pull(pid)
pid_swing2 <- df %>% subset(., !is.na(lrscale) & syear %in% c(2014) & lrscale %in% c(3:7)) %>% pull(pid)
df <- df %>% 
  mutate(swingvoter1_yes = as.numeric(pid %in% pid_swing1), swingvoter2_yes = as.numeric(pid %in% pid_swing2),
         swingvoter1_no = as.numeric(!pid %in% pid_swing1), swingvoter2_no = as.numeric(!pid %in% pid_swing2))

mediatorvec <- c('swingvoter1_yes', 'swingvoter2_yes', 'swingvoter1_no', 'swingvoter2_no')
slope_temp <- "separate"
order_temp <- 1

reg <- list()
dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.))
reg_baseline <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp)
for(m in 1:length(mediatorvec)){
  mediator_temp <- mediatorvec[m]
  dat_temp1 <- rdd_data(df$y, df$x, cutpoint = cutoff) %>% subset(., complete.cases(.))
  df1 <- df %>% subset(., get(mediator_temp) == 1)
  dat_temp1 <- rdd_data(df1$y, df1$x, covar = df1[,contVar], cutpoint = cutoff) %>% subset(., complete.cases(.))
  reg[[m]] <- dat_temp1 %>% rdd_reg_lm(slope = slope_temp, bw = bw, order = order_temp)
  reg[[m]]$pval <- t.test2(reg_model1 = reg_baseline, reg_model2 = reg[[m]]) %>% .[['p.value']] %>% round(., 4)
  reg[[m]]$baseline <- mediator_temp
}
fileName <- paste0('results/APPENDIX_B/B14.html')
names(reg) <- seq(1, length(reg), 1) %>% paste0('(', ., ')')
info <- data.frame() %>% 
  rbind(c('p-value', unlist(lapply(reg, function(x) x$pval)) %>% as.character())) %>%
  rbind(c('sub-sample', unlist(lapply(reg, function(x) x$baseline)) %>% as.character()))
modelsummary(reg, vcov = 'HC1', exponentiate = F, 
             estimate = "{estimate}{stars}", fmt = 3, stars = stars, add_rows = info,
             coef_map = coef_map, gof_map = gof_map, 
             notes = list(note2, note1), output = fileName)

